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iBfiARY 



ABSTRACT 

A finite element formulation for solving axisymmetric 
transient heat conduction and thermal stress problems is 
developed in this thesis. The governing equations of 
uncoupled, linear, isotropic thermoelasticity are discretized 
using quadratic isoparametric elements. A FORTRAN IV program, 
using double precision arithmetic, is presented. Compact 
storage techniques for banded symmetric matrices are used. 

Comparisons between exact and computer solutions demon- 
strate close agreement for a number of test problems. De- 
tailed instructions for using the program are included. 
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LIST OF SYMBOLS 



Note: A single underline is used to denote a column vector 

and a double underline denotes a rectangular matrix. 
The symbols used in computer program are described in 
the beginning of Appendix B. 

a Entrance fluid cross-sectional area 

A Linear combination of C and Y matrices 

b A constant vector 

B Standard rectangular strain- displacement matrix 

C Thermal capacitance matrix 

c Specific heat 

£ Standard elasticity matrix 

E Young's modulus of elasticity 

e Superscript designating element contribution 

F Load vector 

F,p Element thermal load vector 

0 

Fp Element pressure load vector 

F c Element centrifugal load vector 

£ Linear combination of £ and Y matrices 

h Surface heat transfer coefficient 

^ Identity matrix 

J Jacobian coordinate transformation matrix 

K System stiffness matrix 

K Element stiffness matrix 

k Thermal conductivity 

L Thickness of slab 
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I Arc length along the side of quadrilateral 

Shape function 

m Total number of nodes 

n Outward normal or number of nodes per element 

P Pressure 

R Radial coordinate 

S Surface area 

T Nodal temperature vector 

T Temperature or, when used as a superscript, 

transpose of a matrix 

T aV g Average temperature 

T^ Fluid temperature (used in one-dimensional example) 

T 

U A vector defined as <1 1 1 0> 

u Radial displacement 

v Right-hand side vector in conduction equation 

V Volume 

W Work done by loads 

W Modal matrix 

w Axial displacement 

w Eigenvector 

Y,Y Thermal admittance matrix 

y Dependent variable 

* 

y^j Element ij of the matrix Y* 

t ] Matrix representation 

< > Row vector 

V Gradient operator 

a Coefficient of thermal expansion 

8 Constant coefficient vector 
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Nodal displacement vector 
6_ e Element displacement vector 

£ Strain vector 

y Eigenvalue of one -dimensional transient temperature 

solution 

e Thermal strain vector 

— o 

£ e Element strain vector 

T) Local element coordinate 

v Poisson's ratio 

X Eigenvalue 

p Material density 

o Stress vector 

E, Local element coordinate 

x Time 

At Step size of numerical time integration 

t R £ Shearing stress component 

0 Fluid temperature, or angle 

A Spectral matrix 

ft Speed of rotation 
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I. INTRODUCTION 



Thermal stresses have become increasingly important in 
engineering practice during recent years. In power genera- 
tion higher cycle temperatures and use of nuclear fission 
are largely responsible for this trend. This thesis describes 
a computer program for finding temperatures and stresses in 
bodies having axisymmetric geometry and loading. The govern- 
ing equations are those of isotropic, uncoupled, quasi-static, 
linear thermoelasticity. They are discretized by using the 
finite element method. A FORTRAN IV computer program using 
double-precision arithmetic has been written to solve 
problems of the following kinds. 

A. TEMPERATURE PROBLEMS 

The transient temperature vector, evaluated at the nodal 
points , may be obtained for an axisymmetric body with the 
combination of insulated, convection, or constant temperature 
boundary conditions. For the convection thermal boundary 
condition, however, we may have fluid flowing with entry 
temperature prescribed as a linear function of time (RAMP). 

The program can handle up to 15 different ramps, each having 
a different flow velocity, and with discontinuities between 
successive ramps. 

B. STRESS PROBLEMS 

The program will generate load vectors for pressure load- 
ing, centrifugal loading, and axial force. Provision is 
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made for direct input of one additional load vector. Stresses 
may be found for any combination of these loadings. 

C. THERMAL STRESS PROBLEMS 

Thermal stresses may be found for as many as 20 different 
temperature vectors which may be output of the temperature 
problem or direct input. In short, in this part any combina- 
tion of the temperature and stress problems may be used. 
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II. FINITE ELEMENT FORMULATION OF HEAT CONDUCTION 
IN AXI SYMMETRIC BODIES 



A. METHOD OF FORMULATION 

For bodies of revolution under axisymmetric loading the 
mathematical problems presented are two-dimensional. The 
governing equation for non-steady heat conduction is 



where k is the thermal conductivity, p the density, c the 
specific heat, T the temperature and V the gradient operator. 
The superior dot denotes a time derivative. 

Applying Galerkin's principle [1] gives 



where the integral is over the volume V of the conducting 
body and N^ is a "shape function" used in representing the 
temperature distribution. If S is the surface of the body 
and n the outward normal to surface, then using Gauss' 
theorem we can write 



V • kVT = pcT, 



( 1 ) 



/ N.V-kVT dV = / pc N. T dV, 
V 1 V 1 



( 2 ) 



/ V * (N . kVT) dV = / N. k^- dS. 
V 1 S 1 9n 



(3) 



Since 



/ V* (N- kVT) d V = / (VN.) • (kVT)dv 
V 1 V 1 



+ / N. V- (kVT)dV, 
V 1 



( 4 ) 



we can combine Eqs . 2, 3 and 4 to get 
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( 5 ) 



/ pcN . T dV + / (VN . ) • (kVT) dV = / N . k-jff dS . 

V 1 V 1 S 1 n 

Each node of the solid region has a separate discretized 
linear equation calculated from Eq . 5 using the appropriate 
shape function N^. Thus each of the volume integrals on the 
left hand side of Eq . 5 yields a square coefficient matrix 
in the assembled set of equations. Calculation of these 
matrices is a standard process. Details are given by 
Zienkiewicz [ 1 ] . 

The discretized set of equations takes the form 

C T + Y T = v , (6) l/ 

where there is a term by term correspondence with Eq . 5. 

The real symmetric matrices £ and Y represent, respectively, 
the thermal capacitance and thermal admittance. The elements 
of the vector T are nodal temperatures. The vector v, dis- 
cussed in the following section, depends upon the thermal 
boundary conditions. 

In the present development piecewise constant material 
properties have been assumed, i.e., k, p and c are constant 
within each element, but may vary from element to element. 
Also two-dimensional isoparametric elements are used. Appli- 
cable equations are summarized in Appendix A. 



/ 

In this text a double underline denotes a rectangular 
matrix and single underline denotes a column vector. 
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B. THERMAL BOUNDARY CONDITIONS 



Thermal boundary conditions affect only those scalar 
equations of Eq. 6 which correspond to boundary nodes. 
Accordingly, the vector v is sparse. Also, in a single 
problem it is common to have different thermal boundary con- 
ditions on individual portions of the boundary. In what 
follows the subvectors of v (distinguished by individual 
superscripts) which correspond to separate boundary condi- 
tions are treated individually. 

1. Insulated 



ponding to this boundary condition is zero, since = 0 . 
2. Convection 



of the solid and fluid. If the fluid temperature is 0 and 
the heat transfer coefficient is h, then equating heat con- 
ducted away from the surface to the efflux from the solid 
[2] gives 



where the subscript S means that the derivative is evaluated 
at the surface. 



It is clear that for insulated boundary conditions 
the subvector v^^ of the right hand side of Eq. 6 corres- 



3T 



The heat transfer mechanism occurs in the interface 




(7) 



For constant h: 




dS = h / N. (0-T) dS. 
S 1 



( 8 ) 
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In what follows the fluid temperature 0 is taken to be a 
specified function of position and time. For purposes of 
discretization, the fluid temperature is specified at a 
discrete number of fluid "nodes." If 0^ represents the 
fluid temperature at fluid node j , then the fluid temperature 
along the boundary may be represented by 



0 = Z N.0. , 
J 1 ’ 



( 9 ) 



where the Nj 
used for the 



are one-dimensional forms of the shape 
solid. Substitution in Eq. 8 gives 



functions 



i N i k i dS ' * rljfj-V 

^ J 



CIO) 



where the summation extends over the surface nodes and the 



coefficients y^^ are given by 



>Vj ■ h ' N i N j ds 



( 11 ) 



Assembling the contributions from Eq . 10, the subvector 

( 2 ) 

v for the convection boundary condition may be written 



v 



( 2 ) 




( 12 ) 



The contributions - E y..T. 

J 



from Eq. 10 are included by 



augmenting the matrix Y (see Eq. 13 below). 

3 . Constant Temperature 

If 0 represents the constant temperature desired at 
the wetted surface, then we can use the convection boundary 
condition and replace h by a big number (say 10 ). Since h 

is very large, then for thermal equilibrium the temperature 
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T at the surface will be forced to equal 0. So the subvector 

f 31 

v v 1 for the portion corresponding to the constant temperature 
boundary condition can be obtained from Eq . 12. 

Upon the application of these boundary conditions in 
a single problem, the right-hand side vector v will be com- 
bined from the corresponding subvectors and the finite ele- 
ment discretized equation becomes 

C T + Y + T = v, (13) 

where Y + = Y + Y* . 

C. EXACT TIME SOLUTION WITH SPATIAL DISCRETIZATION 

We consider only the solution of Eq. 13 for v = constant 
with T = a at time x = 0. Let T g be a particular solution 
(steady state) with T s = 0 so that 

l s = (iV 1 X* (14) 

For the homogeneous equation 

C T + Y + T = 0 (15) 

the assumption T = w exp(-Xx), where w is a vector and X is 
a scalar, yields the form 

Y + w=X s Cw. (16) 

It is apparent that Eq . 16 defines an eigenvalue problem. 

Let A be the spectral matrix and W the modal matrix with 
normalization according to 

W T C W = I , (17) 
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where ^ is the identity matrix of the same order as C. Now 
let 

T = W exp(-Ax) b (18) 

where b is a constant vector. Substituting this in Eq . 15 
gives 

Y + W exp (-Ax) b = C W A exp (-A t) b. (19) 

Now Eq . 19 is satisfied for all b if 

W T Y + W = A , (19') 

and this is guaranteed to be satisfied since A and W are 
spectral and modal matrices for the eigenvalue problem of 
Eq. 16 with W normalized according to Eq . 17. 

Returning to the original problem (Eq. 13), the com- 
plete solution may be written as 

T = W (B_ + exp(- : AT) b) (20) 

where 

! s = W 6 , (21) 

and B is a constant vector. Now £ may be found (using 
Eq. 14) to be 

$ = A" 1 W T v . (22) 

Substituting this result into Eq . 20 and using the initial 
condition gives the result 

b=W^£-A^ W T v . (23) 

The general solution of Eq . 13 may thus be written 
as 
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T = W(I - exp (- Ax)) A -1 W T v + W exp (-Ax)W" 1 a. (24) 

For purposes of the present program non- zero compo- 
nents of vector v are to be specified as piecewise linear 
functions of time. During each segment of time history of v 
an analytical solution of Eq. 13 is possible in the form of 
a particular solution plus a complementary solution such as 
Eq. 20. At each node the corresponding time variation of 
temperature will consist of a linear part contributed by the 
particular solution and a sum of n terms representing the 
complementary part. Each of these n terms decays exponen- 
tially with a separate time constant. In principle it is a 
straightforward process to find each particular solution and 
accompanying complementary solution. 

Contemplated problems may typically have from 10 to 
40 segments required to represent the piecewise linear 
variation of v. The number of body nodes n will be of the 
order of 100 or more. In view of the number of particular 
solutions required, each accompanied by an individual comple- 
mentary solution of the form given by Eq . 20, it was concluded 
that a numerical solution of Eq. 13 would be considerably 
more economical than an analytic one such as that given by 
Eq . 24 . 

D. TIME INTEGRATION 

In this section the relative merits of the Runge-Kutta 
and trapezoidal methods of time integration are discussed. 
Since either of these methods will give an exact result if 
the solution is a linear function of time, investigation 
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is confined to performance on a single scalar equation 

y + Xy = 0 (25) 

whose solution y = y Q exp(-Xr) is of the same form as the 
components of the complementary solution (Eq. 20 ). 

1 . Runge-Kutta Method 

A method introduced by Runge and subsequently elabo- 
rated by Heun and Kutta [3] is widely used for the numerical 
solution of first order ordinary differential equations. 

This algorithm prescribes a sequence of calculations for 
determining the ordinate y. , at time t • = t-+At in terms 

6 1 i+ 1 l+l l 

of y^ and values of y at intermediate and end points of the 

interval At. The fourth-order form, which requires four 

evaluations of y, gives for Eq. 25 the result 

lili . 1 - XiT . iMlli - . (26) 

The right-hand side of Eq. 26 represents the first five terms 
of the Taylor expansion of the exact solution 
^i+l^i = exp(-XAx)) so we may conclude that the relative 
error in each time step is less than modulus of the next 
term: (XAt) ^/5 ! . 

In addition to providing the apparent prospect for 
high precision indicated by this error bound, the Runge-Kutta 
method also permits changes of time increment during the 
integration process without requiring additional computa- 
tionally expensive matrix decompositions. The attractiveness 
of these two features dictated a thorough exploration of the 
potential of this method for the present application. The 
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is confined to performance on a single scalar equation 

y + Xy = 0 (25) 

whose solution y = y Q exp(-Xt) is of the same form as the 
components of the complementary solution (Eq. 20 ). 

1 . Runge-Kutta Method 

A method introduced by Runge and subsequently elabo- 
rated by Heun and Kutta [3] is widely used for the numerical 
solution of first order ordinary differential equations. 

This algorithm prescribes a sequence of calculations for 
determining the ordinate y^ + ^ at time x^ + ^ = t^+At in terms 
of y^ and values of y at intermediate and end points of the 
interval Ax. The fourth-order form, which requires four 
evaluations of y, gives for Eq . 25 the result 

liii = 1 - X4 t * * OArP _ (26) 

The right-hand side of Eq . 26 represents the first five terms 
of the Taylor expansion of the exact solution 
(y^ + ^/y^ = exp(-XAx)) so we may conclude that the relative 
error in each time step is less than modulus of the next 
term: (XAx) ^ / 5 I . 

In addition to providing the apparent prospect for 
high precision indicated by this error bound, the Runge-Kutta 
method also permits changes of time increment during the 
integration process without requiring additional computa- 
tionally expensive matrix decompositions. The attractiveness 
of these two features dictated a thorough exploration of the 
potential of this method for the present application. The 
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disqualifying defect which emerged after studying a number 
of examples is readily appreciated from examination of 
Table I. For values of AAx less than 0.5 it is apparent 
that Runge-Kutta scheme affords acceptable engineering 
accuracy. However, when the method is applied to solution 
of Eq. 13 we must deal with a number of A's equal to the 
number of nodes (see Eq. 20 ). This number may be greater 
than 100 and the ratio of the largest A to the smallest may 
easily exceed 1000. Although the solution is dominated by 
the contributions of the eigenvectors corresponding to the 
smaller As, it is clear that the solution will be unstable 
if the largest AAr exceeds about 2.7. Because an unaccept- 
ably small At is required in typical problems, the Runge- 
Kutta method was rejected. 

2 . Trapezoidal Method 

The trapezoidal method estimates y^ + ^ from the formula 



y i*l ” y i * + y i + l)- 

Substituting for y^ and y^ + ^ from Eq. 
gives 



(27) 

25 and rearranging 



= I 7 vv T • (28) 

y i 2 + AAx v j 

Series expansion of the right-hand side provides an error 
bound (per s tep) : 

(AAx) 3/12 . 

From the point of view of the size of AAx this method has no 
stability limit, but has slow attenuation with alteration 
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TABLE I 



ATTENUATION FACTOR COMPARISON 
Fourth Order Runge-Kutta Algorithm 



AAt 


(Runge - Kutta) 


exp (-AA t) 
(Exact) 


. 0001 


. 9999 


. 9999 


. 001 


.9990 


.9990 


.01 


.9901 


. 9901 


.1 


.9048 


.9048 


.2 


. 8187 


.8187 


.5 


. 6068 


.6065 


1.0 


.3750 


. 3679 


2.0 


. 3333 


.1353 


2.5 


. 6484 


. 0821 


3.0 


1.3750 


.0498 


4.0 


5.0000 


.0183 


8.0 


110.3333 


. 0003 


10.0 


291.0000 


.0000 


20.0 


5514.3333 


.0000 


50.0 


240784.3333 


.0000 


100.0 


4004901.0000 


.0000 
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in sign for large AAt. Table II shows this behavior. 

Since a wide usable range of AAt is essential and 
the stability of trapezoidal integration is guaranteed, 
this method is chosen for the present program. 

Applying the trapezoidal algorithm to Eq. 13 yields 

A T 1+1 = £ T 1 + y (v 1+1 + v 1 ) , (29) 

where 

A = £ + y Y ' f , 

G = £ - y Y + . 

and the superscripts denote evaluation at discrete time 

intervals Ax. If m is the order of the capacitance and 

admittance matrices, £ and Y, then once a certain step size 

3 

Ax is chosen, it requires m /3 operations to perform the 
needed triangular decomposition of A. Thus, for large m, a 
change of step size At becomes costly from the point of view 
of computer time. Accordingly, in the present program only 
one time step size is used throughout each problem. 

Also, for assuring sufficiently rapid attenuation of the 
components corresponding to the large AAt, the following 
correction is utilized. 

3 . Irons' Correction 

Irons proposed a scheme [4] for augmenting the atten- 
uation of the contributions of those eigenvectors for which 
AAt is large. Define 
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TABLE II 



ATTENUATION FACTOR COMPARISON 
Trapezoidal Integration 



AAt 


y i+i /y i 

(Trapezoidal) 


- XAt 
e 

(Exact) 


.0001 


.9999 


.9999 


.001 


.9990 


.9990 


.01 


.9901 


.9901 


.1 


.9048 


.9048 


.2 


. 8182 


. 8187 


.3 


. 7391 


. 7408 


.5 


.6000 


.6065 


1.0 


. 3333 


.3679 


2.0 


.0000 


.1353 


2.5 


-.1111 


. 0821 


3.0 


-.2000 


.0498 


4.0 


-.3333 


.0183 


8.0 


-.6000 


.0003 


10.0 


-.6667 


.0000 


20.0 


-.8182 


.0000 


50.0 


-.9231 


.0000 


100.0 


-.9608 


.0000 
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( 30 ) 




where y^ and Y^+i a re obtained from y. , by trapezoidal 
integration . 

In the program presented in Appendix B, Eq. 30 is used 
after every 10 steps of time integration. Table III shows 
the resulting modifications. 



E. ESTIMATION OF EXTREME EIGENVALUES 

Analytic results for one-dimensional heat conduction 
give, for an eigenvalue, 




where d is the distance between points of extreme tempera- 
ture and zero temperature. 

If we use this to estimate the smallest A in cylindrical 
coordinates, two modifications are recommended. 

1. Assume that the point of zero temperature is in 
the fluid at a distance from the wall equal to k/h, 
where h is the surface heat transfer coefficient. 

2. If there are two approximately orthogonal paths for 
heat flow from the (single) maximum temperature point, 



coefficient has no significant effect. We may continue to 



2 

then replace 1/d in the above formula by 



1 




(32) 



d 



2 



For estimating the largest A, the surface heat transfer 
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TABLE III 



EFFECTS OF USING IRONS' CORRECTION 
AFTER 10 STEPS OF INTEGRATION 



AAt 


exp (- IOAA t) 
Exact 


yio^o 

Trapezoidal 


7 10^ y 0 
Corrected 


0.00010 


0.99900 


0.99900 


0.99900 


0.00020 


0.99800 


0.99800 


0.99800 


0.00100 


0.99005 


0.99005 


0.99005 


0.01000 


0.90484 


0.90484 


0.90486 


0.10000 


0.36788 


0.36757 


0.36849 


0.20000 


0.13534 


0.13443 


0.13579 


0.30000 


0.04979 


0. 04866 


0.04978 


0.50000 


0.00674 


0.00605 


0.00645 


1.00000 


0.00005 


0.00002 


0.00002 


2.00000 


0.00000 


0.00000 


0.00000 


4.00000 


0.00000 


0.00002 


-0.00001 


8.00000 


0.00000 


0.00605 


-0.00040 


10.00000 


0.00000 


0.01734 


-0.00072 


20.00000 


0.00000 


0.13443 


-0.00136 


50.00000 


0.00000 


0.44914 


-0.00072 


100.00000 


0.00000 


0.67028 


-0.00027 
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use the same formula for -hr , but now consider only the 

d^ 

smallest element and take 



^ _ length of smallest side 

min 4 

^ _ length of largest side 

max 4 



(32’) 



A comparison of estimates based on Eq. 31, 32, 32' with 
the exact solution of the eigenvalue problem has been carried 
out for several examples. Based on these comparisons, it is 
believed that these estimates are sufficiently accurate for 
choosing a time step and estimating the time of occurrence 
of the extreme stresses. 
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III. ONE-DIMENSIONAL HEAT CONDUCTION 



For comparison of numerical (time) integration methods, 
studies of one- dimensional heat conduction were made. In 
this section numerical results for trapezoidal time integra- 
tion using Irons' correction are compared with the exact 
transient temperature solution. 

Consider a flat slab of thickness L with conductivity k, 
density p, specific heat c, zero initial temperature, one 
face insulated and the other in contact with fluid at tem- 
perature T£ (Fig. 1). The surface heat transfer is h. The 
exact transient heat conduction solution is available [5] as 



T = T, 



00 

1 - 2 Z 
i=l 



Siny^L Cosy^x 
y-L+Siny-L Cosy-L 

l l l 



2 

- y i 



pc 



(33) 




Figure 1. Slab with One Face Insulated and Another 
in Contact with Fluid. 
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where k/pc is the diffusivity and the y^ are the solutions 
of the transcendental equation 



Tan yL = 



yL 



(34) 



For the finite element comparisons we subdivide the 
distance L into m one-dimensional 3-noded elements and use 
the corresponding isoparametric shape functions. (The work- 
ing equations, shape functions and the element capacitance 
and admittance matrices are given in Appendix A.) 

In the course of this investigation separate computer 
programs were written to evaluate nodal transient tempera- 
tures using the following methods: 

(a) exact transient temperature solution, Eq. 33; 

(b) exact time solution with spatial discretization 
(section II. C, Eq. 24); 

(c) Runge-Kutta time integration (section II. D, 
part a) ; 

(d) trapezoidal time integration (section II. D, 
part b) ; 

(e) trapezoidal time integration with Irons' 
correction (section II. D, part c) . 

For an initial step change of fluid temperature from 
zero to 1, transient temperatures have been found. For these 
comparisons the parameters (in consistent units) were taken 
to be : 



L= 8, p= 25, k= 8, c = 5 and h = 5 

The distance L was subdivided into m = 3 elements. For the 
present purpose, comparison is confined to the exact solu- 
tion (item (a)), and the finally adopted system (item (e)). 
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In Table IV, temperatures at the two faces (x = 0, x = 8) 
are compared for various times. The trapezoidal integration 
has been performed using the constant time increment unity 
and the Irons' correction is applied after every 10 incre- 
ments . 

It is believed that Table IV demonstrates that the 
numerical integration method gives adequate accuracy for 
engineering applications. 



TABLE IV 

COMPARISON OF ONE-DIMENSIONAL TRANSIENT TEMPERATURES 



Method 


X 


Time= 0 


T ime= 10 


Time=20 


T ime= 30 


Time=50 


T i me = 7 0 


Exact 


0 . 


.00000 


.00000 


. 00000 


.00002 


. 00095 


.00560 


Trapezoidal 
with I rons ' ^ 


0 . 


. 00000 


.00028 


.00007 


.00066 


.00207 


.00619 


Exact 


8. 


.00000 


. 38431 


.47684 


. 53284 


.60264 


.64685 


Trapezoidal 
with I rons ' J 


8. 


. 00000 


. 38720 


.47716 


.53262 


.60256 


.64686 



^Irons' correction is used after every 10 steps of trapezoidal 
integration . 
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IV. STRESS PROBLEM 



For bodies of revolution deformed symmetrically under 
axisymmetric loading, the stress distribution is two- 
dimensional. Since deformation is symmetric about the axis 
of revolution, cylindrical coordinates (R,Z,0) are used. 

It follows that the stress components are independent of the 
angle 0 and all derivatives with respect to 0 are zero. 

Also the components of shearing stress and t^q vanish 

on account of the symmetry. But since any radial displace- 
ment induces a strain in the circumferential direction, 
this non-zero component of strain and the three in-plane 
components (e^, e^, y^) , complete the state of strain at a 
point in any axisymmetric situation. Hence the state of 
stress for an axisymmetric body under axisymmetric loading 
is given by 



- = <a Z a R a 0 t RZ > ' 

In this chapter the stiffness matrix of an axisymmetric 
body and the thermal, pressure, and centrifugal load vectors 
are formulated and, finally, evaluation of stresses at a 
point is discussed. The treatment closely follows that of 
Zienkiewicz [1] and this reference should be consulted for 
further details. 
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A. STIFFNESS MATRIX 



The elements used are bodies of revolution (about the 
Z axis). For analysis it is sufficient to describe the 
cross-section in the R,Z plane. In Fig. 2 such an element 
and the local £,ri coordinates are shown. 




Z 



Figure 2. Quadrilateral Element Representation. 



If u and w are the displacement components at a point in the 
directions of R and Z respectively, then these displacement 
components may be defined in terms of the nodal displacements 
by the appropriate isoparametric shape functions as 



8 8 

u = E N-u-, w = £ N.w., (36) 

i=l 1 1 i=l 1 1 



where N^, a function of £ and n , is the shape function for 
element node i and u^,w^, are the nodal displacement compo- 
nents. The strain-displacement relations [6] can now be used 
to obtain the components of the strain vector. Thus 
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(37) 



e = B 6 



where B is the standard rectangular strain-displacement 
matrix of any finite element formulation, a function of 
the local coordinates Z, and q, and &_ is the vector of nodal 
displacement. (See Appendix A, part 3, where the applicable 
formulas and useful equations are summarized) . 

If the elasticity matrix for an isotropic material is 
D, then the stress vector £ at a point is given by 



Now, by evaluation of the total strain energy in the element, 
the element stiffness matrix can readily be obtained as 



where the integration extends over the volume of the 
element . 

In the present program the upper triangle of each ele- 
ment stiffness matrix is evaluated by numerical integration 
using four Gauss points within the range of E, and of q [1]. 
The element contribution is immediately placed in the system 
stiffness matrix, which is stored in banded form. 

B. THERMAL LOAD VECTOR 

If we denote T as the difference between local tempera- 
ture and reference temperature, then the thermal strain 
is given as 



a = D e 



(38) 



K e = / B T D B dV 



(39) 




( 40 ) 
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I 



where 



U = <1 1 1 0> 



T 



and a is the coefficient of thermal expansion. 

0 

The thermal load vector F^, is given by Zienkiewicz [1] 



as 



It = / = T =2 -o dv • 



(41) 



From the point of view of the numerical evaluation it is 



where E is the modulus of elasticity and v is Poisson's 
ratio. Thus 



where T = E N-T- and n is the number of nodes in each 
i = l 1 1 

element . 

In the attached computer program in Appendix B the advan- 
tage of the simplicity of Eq . 42 has been utilized. Also, 

since B has some zero components, in the process of multi- 
T 

plication of B U, simply the addition of the appropriate 

T 

non-zero components of each column of the B has been per- 
formed. Finally, Eq. 43 has been integrated with four Gauss 
points . 



interesting to note, however, that the product D in 
Eq. 41 will reduce to 




(42) 



-T = i“fv / T = T I dV ’ 



(43) 



n 
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C. PRESSURE LOAD VECTOR 



Consider a quadrilateral element as in Fig. 3 on the 
boundary of the axisymmetric cross-section where constant 
pressure P is applied. The infinitesimal force dF due to 
the normal pressure acting on the inner infinitesimal cir- 
cumferential surface dS is 

dF = PdS = 2 ir RP dl , (44) 

where dii is the infinitesimal length along the side of 
quadrilateral . 




p- Z,w 

Figure 3. Boundary Element Under Pressure. 



Let F R and F^ be the components of the pressure force in the 
R and Z directions respectively, then 

dF R = dF Cos 0, dF z = -dF Sin 6 , (45) 

where Sin 0 = ^ and Cos 0 = . Therefore 

dF R = (2 tt RP) dZ , dF z = -(2irRP)dR. (46) 
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Since the total work W done by the normal force F is equal 
to the sum of the work done in R and Z directions, 

W = / u dF^ + /w dF^ 

= 2ttP (/ R u dZ - /R w dR). ' (47) 

Now, by using the appropriate shape functions, each compo- 
nent of Eq . 47 may be defined in terms of the nodal values. 
Since 

T 

W = (6 e ) F®, (48) 

where F® is the element pressure load vector 
vector contributed by node j is explicitly given as 







r 3N. -i 


f r 


e 

1 

= 2nP / (EN.R.) 


E H Z i 


f 7 


-1 1 1 


3N . 


zj 




L- 1 H R J 



N j d£ 



(49) 



where in Eq. 49 is the appropriate shape function for 
node j . 



D. CENTRIFUGAL LOAD VECTOR 

Refer again to Fig. 2 and assume that the body is rotat- 
ing about the Z axis with angular speed ft. Then the centri- 
fugal force per unit volume at a point distant R from the Z 
axis will be pft R, where p is the density of the material. 
The work done in this case is 

W = / pft 2 R u dV . (50) 

V 
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For constant pfi 2 the corresponding element centrifugal load 
vector is readily obtained by evaluation of the components 

of the integral in Eq. SO in terms of the nodal variables, 

i . e . , 



F^ = 2mp(2 2 (IN.R.) 2 det j Nj dl-dn, (51) 



where det J is the determinant of the Jacobian coordinate 
transformation matrix (see Appendix A). 

E. STRUCTURAL BOUNDARY CONDITIONS 

The structural boundary conditions implemented in the 
program are: 

(a) one or more nodes prevented from moving axially; 

(b) one or more nodes prevented from moving radially; 

Cc) the right-hand end cross - i nnoi 

end boundary condition. S the P lane ‘ 

computer program presented in Appendix B has the 
capability of handling any combination of the above men- 
tioned structural boundary conditions. For boundary condi- 
tions of the types (a) and (b) , simply multiplying the 
corresponding diagonal component of the stiffness matrix by 
10 gives zero displacement [8] (for practical purposes). 

For the boundary condition of type (c) both ends are initially 
axially for all solutions. An additional solution is 
obtained for unit axial displacement of one end. The axial 
force is evaluated for each solution. Superposition is per- 
formed by adding the displacement vectors for the given 
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loadings (thermal plus mechanical) , plus an appropriate 
fraction of the vector found for unit axial displacement. 

This fraction is chosen so that the resultant axial load has 
the specified value. 

F. SYSTEM EQUATION SOLVER 

Once the desired structural boundary conditions are 
applied, then the problem is to find the nodal displacement 
vector 6_, corresponding to a given number of load vectors. 

We have 

£ 6 = F , (52) 

where K is the system stiffness matrix in banded form and F 
is a load vector. In the present computer program a single 
Cholesky decomposition is performed on K. Then, by a process 
of forward and back substitution, each load vector is 
replaced by the corresponding displacement vector. 

G. PRINCIPLE OF SUPERPOSITION 

Upon the evaluation of the displacement vectors due to 
the various types of loading, the principle of superposition 
can be applied on the displacement vectors. On each thermal 
displacement vector the displacement due to any other type 
of loading is superimposed and, as the result, the number of 
displacement vectors is reduced to the number of thermal load 
vectors . 

H. STRESS EVALUATION 

From the system displacement vector <$_, the displacement 
vector of each element <5_ e may be obtained easily. Then the 
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0 

corresponding element strain vector c_ at any point can be 
found from 

£ 5 = B 6_ e . (5 3) 

Finally, the corresponding element stress vector o e is 
obtained by 

£. 6 = D (e e - e®) . (54) 

Since normally the stresses on the inner and outer surfaces 
of axisymmetric bodies are desired, in the computer program 
presented in the Appendix B provision has been made to cal- 
culate the stresses at the two Gauss points corresponding 
to ( = + — — on the inner and outer boundaries of each 

rr 

element (where n = ±1) • 

Upon the evaluation of the stresses at each point the 
mean stress and the octahedral shearing stress [7] are 
calculated. The program gives as output the extreme values 
of these stresses, the R and Z coordinates of the corres- 
ponding points, and the times of occurrence. 
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V. ONE- DIMENSIONAL TRANSIENT STRESSES 



In this section a one-dimensional comparison of stresses 
is made between exact and finite element results. The tran- 
sient temperature problem is the one previously described in 
Section III. 

If the slab edges are free to translate in the plane of 
the slab, but are prevented from rotating, the exact solu- 
tion for thermal stress [9] is 



y 



= a 



Ea 
1 - v 



(T 



avg 



T) 



(55) 



where T is the 

T is 
avg 

If we choose 
units), then the 
tion is 



local temperature (Eq. 33), and 
the average temperature in the slab 

E = 2, a = .50, and v = 0 (all in consistent 
maximum stress obtained from the exact solu- 



o 

max 



-.477786 



and it occurs at x = 8, t = 73. 

Using the finite element technique with trapezoidal time 
integration and Irons' correction every 10 steps, the maximum 
stress is found to be 



cj Y = -.4 77778 
max 



and it also occurs at x 
It is observed that 



= 8, x = 73, .as before, 
the method chosen gives 



excellent 



results . 
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VI. TEST PROBLEMS 



Program integrity and accuracy have been verified by 
solving a number of test problems. Since stresses, whose 
evaluation depends upon derivatives of displacements, are 
known to be less accurate than temperatures, comparisons 
with exact results are confined to stresses. Individual 
problems are described below. 

A. THICK CYLINDER 

Consider a thick cylinder with inside radius 30 inches 
and outside radius 50 inches and the following material 
properties . 

Modulus of elasticity E - 

Poisson's ratio v = 

Coefficient of thermal expansion a = 

Thermal conductivity k = 

Density p = 

Specific heat c = 

An arbitrary length of 25 inches has been 
cylinder and it has been subdivided into two different 
element representations as in Figs. 4 and 5. The plane-end 
boundary condition with zero axial force is used. The 
stresses for various types of loading are compared with the 
corresponding exact analytic solutions as described below. 



28.9 x 10 u Psi 
.28 

7.22 x 10" 6 1/F° 
Btu 



28. 

489. 

.111 



hr. ft°F 
Lbm/ft 3 
Btu 

Lbm . ° F 



selected for the 
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Figure 4. Two Radial Elements Representation 
of Thick Cylinder. 




Figure 5. Five Radial Element Representation 
of Thick Cylinder. 

1 . Thermal Loading 

For a linear variation of the temperature T = 20R 

the stresses obtained by the finite element method for the 

above representations are compared with the exact analytic 

solution in Fig. 6. The for this problem clearly is 

- 9 

zero and the one obtained by the program was 10 . The 

accuracy of the other results is clearly satisfactory. 

2 . Pressure Loading 

A uniform pressure of 1000 psi. acts on the inner 
surface of the cylinder. Again, is zero and the program 

gives 10 In Fig. 7 the other stresses induced by this 

uniform pressure are compared x^ith the exact solutions. 

Here also the accuracy of the results, even with two radial 
elements, is adequate. 
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3 . Centrifugal Loading 

The speed of rotation has been assumed to be 500 
revolutions per minute. The stresses obtained by the pro- 
gram are compared with the analytic solution in Fig. 8. In 
this case also the results obtained by the program, even 
with only two radial elements, are very close to the exact 
solution . 

B. HOLLOW SPHERE 

We consider a hollow sphere with the same material 
properties as in the thick cylinder with inside spherical 
radius 30 inches and outside 50 inches. The loading is 
thermal with T = 20 ’(spherical radius). Symmetry permits 
using only half of the sphere for the computer analysis. 

The elements representation is given in Fig. 9. Since the 
program gives stresses in cylindrical coordinates, these 
have been transformed to the spherical coordinates for 
comparison with the exact solution in Fig. 10. The accuracy 
of the results is noteworthy. 

C. THERMAL STRESSES IN NOZZLE 

This problem concerns thermal stresses near the inter- 
section of a cylindrical pipe and the spherical vessel. 

Fig. 11 gives the cross-section of the structure which is 
to be analysed. The material properties are: 



Modulus 


of elasticity 


E = 29.3 x 10 6 


Ps i 


Poisson ' 


's ratio 


v = .30 


- 


Thermal 


expansion coefficient 


a = 7.6 x 10" 6 


1/F 



46 



KSI. CS" 

4.5 

4 . • 
3 . 5 - 
3 . • 

2.5 

2 . 

45 

/. 

.5 

O — 

-.5 




O 5 elements ■ 



FIG. 8 
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Figure 9. Hollow Semi-sphere 32 Elements Representation. 
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Thermal Stresses in Hollow Sphere 
T = 20 (spherical radius). 
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NOZZLE GEOMETRY 



* <f9J 




K 


o 




£ 


V- 




Uj 
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O 




s 


< 


Uj 


iu 


Uj 


~J 

o 


-J 


U 

o 

Cfc 


ct 

O 


K 

3: 


5 


ct 


o 


o: 


❖ 


£ 

<S) 


*io 

» 

R 


18.078 



O 

lo 

OJ 

c\j 



50 



Thermal conductivity 


k 


= 10.25 


Btu/ft .hr. °F 


Density 


P 


= 530.5 


Lbm/ ft ^ 


Specific heat 


c 


= .128 


Btu/Lbm°F 


loading results from the 


thermal 


transient 


in the fluid 



contained in the nozzle. This fluid is in contact with the 
structure on surface "1" (Fig. 11) and has the entry temper- 
ature time variation as in Fig. 12. The fluid in the sphere, 
which is in contact with the structure on surface "2" (Fig. 
11), has the constant temperature 478°F. At t = 0 structure 
has a uniform temperature of 478°F and is stress free. 
Exterior surface of the structure is insulated. The flow 
velocity past surface 1 is 8.5 ft/sec (inward) and 
there is no flow past surface 2. The surface heat 
transfer coefficients are 1393 and 2910 Btu/hr.ft^°F for 
surfaces "1" and "2" respectively. 

For the structural boundary condition it is assumed 
that the nodes on the left end cross-section of Fig. 13 are 
prevented from any axial displacement. 

The maximum thermal stresses obtained by the program 
occur in element 14 of Fig. 13 as follows: 

Maximum mean stress = 18.81 ksi. at time = 4 seconds; 

Maximum octahedral shearing stress = 11.5 ksi. at time 

= 18 seconds. 

» 

These results appear to be reasonable, but no suitable com- 
parison solution is available. 
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IN NOZZLE. ( SURFACE 7 \) 
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FjC. 12. Fluid Temperature - Time Histories. 
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Element Representation 
in Nozzle . 
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VII. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

A computer program has been developed for the solution 
of axisymmetric transient heat conduction and thermal stress 
problems. The system will accommodate a wide variety of 
geometric arrangements, thermal and structural boundary con- 
ditions, and mechanical loadings. 

Although double-precision arithmetic is employed through- 
out, efficient algorithms for the manipulation and storage 
of large symmetric banded matrices allow in-core solutions 
with modest time requirements. 

The quadratic isoparametric elements used allow accurate 
representation of curvilinear boundaries and the stress 
field. Examples presented show that a small number of ele- 
ments is generally sufficient to determine stresses with 
good precision. 

B. RECOMMENDATIONS 

Incorporation of several additional features would 
significantly increase the program capabilities. The follow- 
ing extensions are recommended. 

1. Material thermal and elastic properties have been 
assumed constant within each element. Since such properties 
are generally temperature dependent, provisions should be 
made for periodically "updating" them during both temperature 
and stress solutions. 
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2. The surface heat transfer coefficient, taken as 
constant in the program, is a function of temperature and 
flow velocity. Provision should be made to include these 
effects . 

3. The program presently starts every temperature 
solution with constant initial solid temperature. Provision 
for an externally specified initial temperature vector should 
be included. 

4. The large quantity of temperature and stress results 
generated by the program is currently presented as digital 
printout. Graphical output in the form of two-dimensional 
contour plots of temperatures and stresses should be 
provided. 
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APPENDIX A 



APPLICABLE FORMULAS AND EQUATIONS 



PART 1 

(a): Two dimensional 8-noded (parabolic) 

isoparametric shape functions. 

Corner nodes : 



N. = ']■ (1 + £ ) (1 + n ) (£ + n - 1) 

i 4 o o o o 



Mid nodes 



C i = °’ N i = \ C 1 ' ? 2 )(! + no- 



where 



r>i = 0, N. = f (1 + 5 0 )C1 



2 -. 
n ) 



» 



£ £ . , n = n n • 
^ ^1 ’ o 1 1 



(b) : Temperature at a point in terms of the 

nodal temperatures. 



8 

T = 2 N.T. 

i-1 1 1 



(c) : Coordinates at a point in terms of the 

nodal coordinates. 



8 

R = 2 N.R. 

i-1 1 1 



8 

Z = 2 N • Z - 

i-1 1 1 



(d) : The Jacobian coordinate transformation 

matrix . 
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(e) : Element of the finite element capacitance 

matrix . 

c1\ = pc / N i Nj dV 

(f) : Element of the finite element admittance matrix. 

y® . = k f VN . • VN . dV 
7 ij i J 

(g) : Elemental volume. 

dV = 2irR det J d£ dp 



PART 2 

(a): One-dimensional 3-noded (parabolic) iso- 

parametric shape functions. 



End 


nodes 


N. = 

l 


I C 

2 o 


(1 - C Q ) 


Mid 


node 


II 

•H 


(1 - 


s 2 ) 


where, again, 


^o = 







(b) : One-dimensional capacitance and admittance 

matrices . 



r e _ pc£ 
= 30 



4 

2 

•1 



2 -1 
16 2 

2 4 



_k_ 

3£ 



7-8 1 

8 16 -8 
1-8 7 



0 0 
0 0 
0 0 



where £ is the length of the element. 



(c) : The element v vector for the one-dimensional case. 



PART 3 



(a): Strain- displacement relations 



£ - 9W 

Z 8 Z 



3N, 

z Tz w i 



_ 3U _ _ 9N i u 
e R 9R ' 1 9 R u i 



I N • U . 

U 11 



e R Z N i R i 



. . 3N. 3N . 

y s iLy. + = y . — u + z i. 

Y RZ 3Z 3R L 3z 1 3r 



w . 
i 



(b) 



where NR are the same isoparametric shape 
function as in Part 1, (a). 

The B matrix. 



B = 



0 

3N-, 



3R 
N , 



R 

3N, 



3 Z 



3N ] 

"3Z~ 

0 

0 

3N ] 

W 



0 

3N. 

w 

N 2 

R 

3N. 

w 



3 V, 
3Z” 

0 

0 

3N, 

W 



3N 



8 



Z 

0 

0 
3 N. 



8 



3R 



(c) : The element displacement vector, 



(d) 



6 = <Uj Kj u 2 w 2 



u s V 



Elasticity matrix D for an isotropic material 



D = 



E(l-v) 

(l + v)(l-2v) 



v 



1-v 



V 



1-V 
V V 



1-v 1-v 
0 0 
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1-v 
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2(l-v) 



APPENDIX B: COMPUTER PROGRAM 
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APPENDIX C 



USER'S MANUAL 

In this section the procedure for the use of the present 
computer program is described. It is intended that a person 
with minimum familiarity with the details of the finite 
element method and computer programming be able to use this 
program . 

The steps below are in order and the user is advised to 
follow them carefully. 

For finding the thermal stresses in an axisymmetric body 
under axisymmetric loading, the user has the choice of using 
either the Metric or British system of units. The program 
will handle both systems and the necessary conversions are 
made automatically within the program. However, .the units 
used in each system must be consistent and they should be as 
listed below in Table V. 

Now for preparation of the data input for a given axi- 
symmetric geometry with prescribed thermal and structural 
boundary conditions we go through the details with a simple 
example . 

Step 1 : 

Draw the longitudinal cross-section of the body to scale. 
Identify the cylindrical coordinates R and Z, with the 
origin of the Z axis on the most left-hand point of the 
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TABLE V 



UNITS FOR INPUT DATA 



Note: an input card specifies whether 'British or Metric 

data is being used. 



Quantity 


British Unit Syst. 


Metric Syst. 


Coordinates 


in . 


cm 


Time 


sec . 


sec . 


Temperature 


° F 


°C 


Velocity 


f t . /sec . 


m/sec . 


Axial force^ 


lbf 


kg 


Pressure 


lbf/in 2 


kg/cm 2 


Density 


lbm/ ft ^ 


gm/ crn"^ 


Specific heat 


Btu/ lbm° F 


cal /kg° C 


Mod. of elasticity 


lbf/in 2 


kg /mm 2 


Coef. of thermal exp. 


1/ °F 


1/°C 


Thermal conductivity 


Btu/hr . ft . 0 F 


cal/sec .cm. °C 


Heat transfer coef. 


Btu/hr. ft. 2 °F 


? 

cal/ sec. cm. °C 


Load vector 


lbf 


kg 


Rotational speed 


Revolutions/min . 


Revolutions/min . 
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body. The entire cross section must lie in the first quad- 
rant of the coordinate plane. 




Figure 14. Longitudinal Cross Section. 

Step 2 : 

Subdivide the cross section into a maximum of 40 quadri- 
lateral elements. This subdivision is extremely important 
and we should provide smaller elements for the parts of the 
cross section where we expect the largest stresses or tem- 
perature gradient. If the body is made from several mater- 
ials, elements should be chosen so that each element contains 
only one material. Any two adjacent elements must share one 
complete side of the quadrilateral. Number the elements, 
starting from 1, in any arbitrary manner (see Fig. 15). 

Step 3 : 

Since eight-noded elements are used in the program, iden- 
tify these nodal points around the boundary of each element. 
Four nodes will be at the corners and the other four will be 
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at the mid-points of the sides. Number sequentially these 
nodes starting from 1 and increasing in the direction where 
the number of elements is the least. (The numbers thus 
assigned are called global node numbers.) For clarity of 
this step assume a cross-section as in Figure 15 such that 
the maximum number of elements in one direction is less than 
the maximum number of elements in other direction. (In 
Fig. 15 we have maxima of 2 and 3 elements in R and Z direc- 
tions respectively.) Thus we number the nodes beginning in 
the R direction. See Fig. 16. This method will give the 
minimum band width of the stiffness matrix and will save 
computer execution time. 




Step 4 : 

At this point we are ready to prepare the first data 
card. Input quantities are: 
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Figure 16. Element and Node Numbering. 

NET The total number of elements 

NNT Total number of nodes 

NCN Total number of corner nodes 

N0FN Total number of mid-side nodes not lying on the 
straight line joining the corner nodes (number 
of "Off" nodes) 

NMAT Number of different materials 

NPR0B Number of problems to be solved 

For the present we take NPR0B = 1. 

Prepare a single card that reads 

NET, NNT, NCN, N0FN , NMAT, NPR0B 

with the format ( 8 1 5 ) . For our example, if all the elements 
are from the same material this data card reads: 

4 23 1 0 111 

0 0 0 C 0 0 C 0 C 0 0 0 0 0 G 0 0 0 0 0 0 0 0 0 0 0 0 0 0 G 0 0 C 0 0 0 0 0 0 G 0 0 C 8 ? fl 3 0 0 0 0 3 G 8 0 C 0 3 8 C t 

1 1 J ( 5 6 i 8 3 13 11 ■ ! 13 is IS 17 IS i« 23 ?i22 23 il 25 25'27 28 :3 30 3) 52 38 )< JS 35 31 38 53<0<l <2 <1 <«<$<: il 'S SD 3' >2:3 51 55 16 5! ae jt i : . 5! s 
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On this card, or any succeeding one, if an input quantity 
(e.g., N0FN) is zero, the corresponding field may have a 0 
or be left blank, unless otherwise stated. 

Step 5 : 

For each corner node a card should be punched which 
reads: The node number I, the R coordinate and the Z coor- 

dinate of the node. The format is (15, 2F10.5). 

It is not necessary to sort these cards in the order of 
increasing or decreasing corner node numbers, they can be 
put together in any arbitrary order. A typical card is 
illustrated in step 7. There must be as many punched cards 
as the number of corner nodes (NCN) . 

Step 6 : 

In this step we read in the connectivity array, i.e., 
for each element we prepare one card which gives the global 
node numbers and the material identification number. 

For each element, start from the lower left-hand node 
and move in the counterclockwise direction within that 
element and punch the global node numbers in order. The 
format is (915) where the last 15 is the material identifi- 
cation number. 

It is very important that the cards prepared for this 
step be put together in order of elements, i.e., the first 
card for element 1, the second card for element 2, etc. 

For ease of sorting the connectivity cards the element 
number may be punched after column 50. As an example, for 
element 2 of Fig. 16 the connectivity card should read: 
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0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 C 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 I) C 0 C C 0 0 0 0 0 C C 0 r c 8 0 i 
! . ? ! 5 ? ! ? >, '»'J»l<l5«'n«lJa?JMm«?5K?>?J2J101l3?13J43S361?3U}«<l«IJ 44 4S4i;<r««SJSIS?SJS«SSSfS75IS!M6II 



1111111111111111^ 



where 1 identifies the type of material in this element and 



2 in column 55 stands for the element number. 



Step 7 : 

If each side of every element is straight, N0FN = 0 and 
this step is omitted. 

For each mid-side node that is on a curved element edge 
a card is prepared with format (I5,2F10.5) to read the global 
node number and the R and Z coordinates. Here, as in step 5, 
these cards may be put together in any arbitrary order. For 
the case of Fig. 16, there will be only one card to be punched 
and it would read: 



i O 
x ■_* 








0000000000000 0000000 OOOOOOCOOODGOOOOOOOOUOGGfi , jO(!OPOOi' i u8 0GC' ! ' 

! ! I 4 ! ? ? 6 ? [“ " [ ! ? '* '. 5 '» [’ '• » n ? l K K U n « !’ 29 31 II 3? 33 31 35 36 37 3: 33 15 <1 1? 43 « 45 46 1. 19 19 50 5' 5? S3 51 5' St 57 56 55 60 El • 
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, To assure in all cases that the double precision capability 
of the program is not compromised by less precise floating 
point inputs, it is recommended that all such inputs be made 
in D- FORMAT. The F-FORMAT instructions merely allocate 
input card fields. 
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Step 8 : 



The properties of the different materials are specified 
in this step. For each material a card must be prepared 
that gives the modulus of elasticity (E) , coefficient of 
thermal expansion (AL) , Poisson's ratio (P0I), thermal con- 
ductivity (TK) , density (DENS) and specific heat (SHT) . 

I . e . , read : 

E, AL, P0I, TK, DENS, SHT 

with the format (2D12.4, 4F8.3). The first card will be for 
material number 1, the second card for material number 2, 
etc . 

Step 9 : 

For this step a single card must be punched, starting in 
column one, which reads the word BRITISH or METRIC in accor- 
dance with the system of units used. 

This completes the input of geometric and material prop- 
erty data. 

Step 10 : 

For a transient temperature problem only (no stress cal- 
culations) leave a blank card for this step and proceed 
directly to step 16. 

For a stress problem or thermal stress problem, in this 
step we specify the type of problem and the structural 
boundary conditions. 

The following quantities, when pertinent, are to be 
specified. 
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0MEGA 



The speed of rotation, about the Z axis, in 
revolutions per minute 

PRES The external pressure applied to a boundary 
segment 

IEXT The number of known temperature vectors for 

which evaluation of thermal stresses is desired 

L0AD The indication for an additional known load 
vector. If there is one, L0AD = 1, otherwise 
L0AD = 0 

NNLT Total number of nodes fixed in the longitudinal 
direction 

NNRT Total number of nodes fixed in the radial 
direction 

IPLANE An indication for the plane-end boundary condi- 
tion. If desired IPLANE = 1 

If any of these items is not applicable, the correspond- 
ing field is left blank. 

Now we have to punch a single card for this step which 
reads : 

0MEGA, PRES, IEXT, L0AD , NNLT, NNRT, IPLANE 
with format (2F10.4, 515). 

Step 1 1 : 

If there is no pressure loading (PRES=0), omit this 
step. For non-zero pressure, the total number of nodes on 
pressure loading boundary segment (NPNT) and the global 
node numbers of these (pressure) nodes (NPN(I)) must be 
specified here. Only one segment is permitted. 

We read in: 

NPNT, (NPN(I) ,1=1 , NPNT) 
with the format (1215). 
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In our example, Fig. 16, 
body, then the data card for 

*71 4 6 11 14 



if there is pressure inside the 
this step will be: 

19 £1 



0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 D C D 0 0 (1 0 0 5 C 0 0 0 e 0 C 0 0 ( 

1 2 3 4 5 6 / 8 9 10 II 12 13 14 15 16 l? 18 19 20 21 22 23 24 25 26 27 2e 28 30 31 32 33 34 35 36 37 38 33 40 41 42 43 44 45 41.47 48 43 50 51 52 53 54 55 56 57 5F 59 bO 61 6 



Step 12 : 

If there is no additional load vector (L0AD = 0), omit 
this step. 

For L0AD = 1 read in the additional load vector with the 
format (6F10.4). The components of the load vector are 
arranged in the order: R component at node 1, Z component 

at node 1, R component at node 2, etc. I.e., READ 
(F (I) , 1=1 ,NDF) , where NDF (the number of degrees of freedom) 
is equal to two times the total number of nodes (2*NNT) . 

Step 13 : 

In any stress problem, or thermal stress problem, there 
must be at least one node constrained against longitudinal 
motion, i.e., NNLT > 0, so we must specify here which nodes 
are to be constrained against longitudinal motion. These 
global node numbers NNL(I) are read in with the format 
(1015), i.e., 

READ: (NN (I) ,1=1, NNLT) 

In our example, if global nodes 1 and 21 are fixed in the Z 
direction, we have a card that reads: 
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1 "21 






oocooooo 

12 3*5616 



ooooooooooooooooooooooo 

9 10 11 12 13 14 is 16 17 18 19 20 21 22 23 24 25 26 27 28 23 30 31 



0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 

32 33 34 35 36 3/ 33 39 40 <1 42 43 44 43 46 V 48 43 50 51 b? S3 54 55 56 3/ 56 5b 60 61 i 

, „ , , , , , ♦ . , , , i ? i « i • ■> • i i » i 1 ' 1 i 1 1 



Step 14 : 

If there are no nodes fixed in the radial direction 
(NNRT = 0), omit this step. 

For NNRT > 0 we must read in the global numbers (NNR(I)) 
of those nodes which are to be fixed against radial motion. 
The format is (1015), i.e., 

READ: (NNR(I) ,1=1 ,NNRT) 

If node 1 and 21 of our example are also fixed in the R 
direction, then the card for this step would read exactly 
as the one in Step 13. 

Step 15 : 

If there is no plane-end boundary condition (IPLANE = 0) , 
omit this step. 

For the case of end-planes -remain-plane boundary condi- 
tion, i.e., (IPLANE = 1), we have to specify the total number 
of nodes of the right-hand end (NNRE) and the global node 
numbers of this end (NNR(I)). 

We read in: 

NNRE, (NNR(I) ,1=1 , NNRE) 
with the format (1015). 

For the case of Fig. 16, if it is desired to have end- 
pl anes - remain-plane boundary condition then the data for 
this step will be: 
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0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 C C 0 0 C 0 0 0 0 0 0 0 0 0 (III 0 3 0 0 0 0 C 0 0 0 C 0 G ( 

1 2 3 4 s 6 7 8 $ 10 l! I? 13 14 15 16 W 16 13 20 21 22 23 24 2$ 25 27 28 29 30 31 32 33 34 3S 3£ JJ 33 73 <3 41 42 <3 44 4*' 4B 47 48 t'l 50 51 52 53 $4 : 5 55 57 58 5$ kfi 61 6 

Step 16 : 

In this step we begin the specification of thermal 
boundary conditions. Additional details are prescribed in 
Steps 19 through 24. 

Thermal boundary conditions are imposed along discrete 
segments of the boundary. Each such segment must begin and 
end at a corner node of a boundary element. We consider 
separately the imposition of the various kinds of thermal 
boundary conditions. 

(a) Convection 

Since the fluid temperature for convection is deter- 
mined from a prescribed temperature - time history at entry 
section and a specified flow velocity (see Appendix D, 

Section 5), it is necessary that the terms "inside," i.e., 
adjacent to the symmetry axis, and "outside" have their 
ordinary meanings. Thus on the inside surface the convec- 
tion boundary condition may be applied to a single (con- 
tinuous) segment. A similar prescription may be employed 
for the outside portion. The entry section used for flow 
calculations is at the upstream end of corresponding segment. 

(b) Constant temperature 

To specify constant temperature on a portion of the 
boundary, the designator "inside" or "outside" may be used. 
Such a constant temperature portion may consist of several 
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discrete segments. If two different constant temperature 
portions are prescribed, one may be designated "inside" and 
the other "outside." 

(c) Combinations of convection and 

constant temperature 

The convection and constant temperature conditions 
can be used together, but the portion designated "inside" 
must have only one of these conditions prescribed and the 
same restriction applies to the "outside." 

(d) Insulated 

All portions of the boundary not included in the 
segments specifically identified as "inside" and "outside" 
are considered insulated. 

The following items, when pertinent, are to be 
given as specified below: 



TINIT The constant initial body temperature 

Q2=0. For insulated boundary condition outside 

Q2 >0. For convection or constant temperature boundary 
condition outside 

Q3=0. For insulated boundary condition inside 

Q3 >0. For convection or constant temperature 
boundary condition inside 

Q4 <0. If solving stress problem only 

Q4 > 0 . If solving temperature problem only 

Q4=0. If solving thermal stress problem 

AXIALF The axial force in the Z direction. As usual, 

+ (plus) for tension and - (minus) for compres- 
sion. 1 



So, we read: 

TINIT, Q2 , Q3 , Q4 , AXIALF 
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with the format (5F10.4). 



In case of Fig. 16, if the initial solid temperature is 
60° and we have insulated outside and convection boundary 
condition inside with 120 kg axial compressive force, then 
the card for this step will be: 

60. DO 1.0 120. DO 

*. uW k Km. . * -r . I* T- -T. < . J*- V .*■ ... * ^ 

0 0 0 0 0 0 0 0 0 0 C 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 n 0 G t 0 0 e 0 0 0 0 0 0 0 fl 0 0 0 0 [ 

1 2 3 < 5 6 r 9 9 10 II 12 13 14 15 16 11 16 19 20 21 22 2 3 24 2 5 2C 22 26 29 30 31 32 32 30.36 3 2 33 35 <0<1<2 43 44 45 4S 42 48 49 5. 51 52 53 54 55 56 5-' 59 59 60 61 5 



where the blanks left in the columns 31 through 40 indicate 
that we want to solve a thermal stress problem. 

Step 17 : 

No action is taken in this step - we merely choose between 
proceeding to step 18 or jumping to step 26. 

If there are no known temperature vectors for which 
evaluation of thermal stresses is desired, proceed to 
step 18. 

For IEXT>0, i.e., when some known temperature vectors 
are to be entered for thermal stress analysis alone or 
combined with some other loadings, proceed directly to 
step 26. 

Step 18 : 

If (Q4<0), i.e., we are to solve only a stress problem, 
proceed directly to step 27. 
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Step 19: 



If there is an insulated boundary condition outside 

(Q2 = 0), omit the following steps and start from step 22. 

For the case of a convection or constant temperature 

boundary condition outside (Q2 >0.), we have to specify 

the outside heat transfer coefficient (HTC0) (for constant 

2 0 

temperature HTC0 =10 .), the constant outside initial 

fluid temperature (TEMP0) which may be equal to the solid's 
initial temperature, the total number of nodes in contact 
with fluid outside (NCF0T) , and the number of temperature 
ramps for outside flow (NRAMP0) . We read: 

HTC0 , TEMP0 , NCF0T, NRAMP0 

with the format (D16 . 4 , F10 . 4 , 21 5) . See example in step 22. 
Step 20 : 

Here we read in the global node numbers of the nodes 
in contact with outside fluid (NCF0(I)). The sequence of 
these node numbers must be in the direction of flow velocity. 
Read : 

(NCF0 (I ) , 1=1 ,NCF0T) 

with the format (1215) (see example in step 23). 

Step 21 : 

For each ramp of outside flow Ave read in the time when 
the ramp starts BDRY0(I,1), the initial temperature of the 
ramp BDRY0(I,2), the final temperature of the ramp (specify 
only if it differs from the initial temperature of the next 
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ramp) BDRY0(I,3) and the velocity of the fluid for this 
ramp BDRY0(I,4), with the format (4F10.4) so we read: 

( (BDRY0 (I , J) ,J=1 ,4) ,1=1 ,NRAMP0) . 

The number of cards prepared at this step is equal to the 
number of ramps. See example in step 24. 

Step 22 : 

If there is an insulated boundary condition inside 

(Q3=0) , proceed directly to step 25. 

For the case of convection or constant temperature 

boundary condition inside, Q3>0. We have to specify the 

inside heat transfer coefficient (HTCJ) (for constant 

2 0 

inside temperature HTCI=10. ), the constant inside initial 
fluid temperature (TEMPI), which may be equal to the solid's 
initial temperature , total number of nodes in contact with 
fluid inside (NCFIT) and the number of ramps of the inside 
flow (NRAMPI). We read: 

HTCI , TEMPI, NCFIT, NRAMPI 

with the format (D16.4, F10.4,2I5). 

For our example assume the time variation of entry temper- 
ature for the inside flow to be as in Fig. 17. The heat 
transfer coefficient is 175.0 with initial inside fluid 
temperature 70.0. The card for this step is: 

175. BO - 70:D0 7 3 

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ? C 0 0 0 0 0 0 0 0 0 0 0 0 0 o C 0 0 0 0 0 0 G G 0 0 0 0 : 

1 l J 1 5 S 1 C 5 13 !! 12 13 14 15 1£ 1113 :5 20 21 22 23 21 23 22 27 21 23 53 3’ 32 3: 'J S'- i: 32 3! 31 12 <1 1; 13 It 131 «■ IS 1} 30 5! 53 5; '1 35 35 5' St r: «0 6if 
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Further details of this temperature- time history are speci- 
fied in step 24. 

Step 23 : 

Here we read in the global numbers of those nodes which 
are in contact with inside fluid, (NCFI(I)) with the format 
(1215). It is important to read these node numbers in the 
direction of flow velocity. 

For example, if the inside flow of Fig. 16 is from left 
to right, then for this step the card reads: 

I 4 6 11 14 1? 2i 

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 9 0 0 0 0 0 0 0 0 0 G 0 G 0 G 0 0 0 G 0 9 0 0 0 0 G 0 0 1 

1 2 3 4 S t 1 8 S 10 111? 13 14 IS IS 17 ‘-8 19 20 2? 23 2* ?S Zi V 28 23 JO 31 32 33 34 35 36 37 3d 33*40 41 42 33 <4 <5 4o 4/ 48 49 50 SI 52 S3 54 SS 56 5/ 58 59 »8 El 8 
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Step 24 : 



In this step data for inside flow are given. For each 
ramp I we read: the time when the ramp starts BURYI (1,1), 

the initial temperature of ramp I BDRYI(I,2), the final 
temperature of ramp I (specify only if it differs from the 
initial temperature of the next ramp) BURYI (1,3), and the 
velocity of the fluid for this ramp BDRYI(I,4). We read in 

( (BDRYI (I,J) , J=1 , 4) ,1=1 ,NRAMPI) 
with format (4F10.4). 

For example of Fig. 17, let the flow velocity for the 
first ramp be 15 and that for the second and third segment 
be 20. Three cards are needed*. 

ISrt.Di'i 3 in Y. DO 300. DO EO.DO 

^ ^ ^ . . v . - . 

ed 

3 CARD 

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 00 0 0000 0000 
1 2 3 i 5 6 7 6 S 10 11 12 U 14 IS 16 17 1* '9 TC 2'. 22 23 21 23 23 27 23 23 35 ?1 3*3. 34 35 3* 37 38 39 40 41 4? 41 14^4 4* !8 5! 52 ‘ 

■ 100. DO 165. DO £00. DO 20. DO 

> &<» <fe 

2 nd CARD 

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 CDG DO 080 8 0 0 0 8 

$ 26 21 ?e 23 30 31 3? 2* h 3: 37 33 39 40 4‘ 42 <j 41 4i It *' O ' «? SG S 1 l 2 S.< << 55 \t 56 

15. DO 

( . . _ ... - .... 

l Sf CARD 

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 OOfifj COO 0 0 0 0 0 0 0 80 0 0 0 01 

1 2 3 4 5 6 7 6 S 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 2$ 27 22 29 33 31 32 32 34 35 36 37 36 33 40 4) 42 43 44 45 46 </ 48 43 50 51 52 S3 54 55 56 57 Sb 58 60 61 F2 63 1 



Step 25 : 

In this step we must specify the following items: 

DTI The time increment (step size) of trapezoidal 

integration 

TIME! The time when the first calculated temperature 
vector is to be stored for thermal stress 
evaluation 



000 00 00000000 POOOOO 

1 2 3 4 5 6 1 : ! 1C II 12 : 3 14 !S1£ If i‘. 13 <D 21 22 23.1 



0. DO 



50. DO 
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EVERY The interval of time between the successive 
temperature vectors which are to be stored 

IVEC Total number of temperature vectors to be 
stored for thermal stress evaluation 

INTP The indication for printing the calculated 
temperature vectors 

If INTP = 0 there will be no temperature prints 

If INTP = 1 the program will print the tem- 
perature after every step of time integration 

If INTP = 5 the program will print the tem- 
perature after every 5 steps of time inte- 
grations, etc. 

At this step we prepare a single card that reads: 

DTI, TIME1 , EVERY, IVEC, INTP 
with format (3F10 . 4 , 2 1 5) . 

For example, if the step size of trapezoidal time inte- 
gration is 2 sec. and we desire to evaluate thermal stresses 
for 10 different temperature vectors, each 40 seconds apart, 
and starting with the first thermal stress calculation at 
time equal to 60, then the card for this step is: 

£,im 60. DO 40. DO 10 15 



UOOOOC 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0(1 0 5 0 0 0 0 0 0 0 0 0 0 0 Of 

1 2 3 4 5 6 7 8 9 10 11 12 13 14 IS 16 17 lb IS 20 21 22 25 24 ?$ 26 27 28 29 30 31 32 V 34 35 35 37 38 39 40 <1 42 43 44 M U </ 48 49 50 51 52 63 54 55 5S 5? '>5 59 60 S! 6 

1 i 1 1 i i < i i i i i i i i ' 1 i ' 1 ' 1 ’ 1 1 * 1 1 * t««i 

where 15 in columns 34 and 35 indicates that the temperatures 
will be printed every 15 steps of time integration. The 
total time of integration for. this example is 420 seconds, 
since 



TIME1 + EVERY * (IVEC-1) = 420 
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Step 26 : 

If there are no known temperature vectors (IEXT = 0) 
for which thermal stress evaluation is desired, omit this 
step . 

For IEXT > 0 we must specify the nodal temperatures 
which are to be stored for thermal stress evaluation. We 
read : 

( (ST0R (I,J) , J=1 , NNT) , 1 = 1, IEXT) 
with the format (6F10.4). 

This means that we enter all the components of the first 
nodal temperature vector, followed by the components of 
the second nodal temperature vector and so on until IEXT 
such vectors have been entered. 

Step 27 : 

The data cards for this problem are completed now. If 
no more problems are to be solved for the same geometry, 
omit this step. 

If another problem is to be solved for the same geometry 
and spatial discretization, increment by 1 the NPR0B in step 
4 and start the input of the new problem from step 10. 

Step 28 : 

If there is no access to the IBM 360 computer of the 
Naval Postgraduate School, omit the following steps and 
start from step 32. 

For the convenience of the user, the author has put 
a so-called CHECK program and the program presented in 
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Appendix B in the computer system at N.P.S. (Naval Post- 
graduate School). It is advised, however, that the user 
check his input data deck with the CHECK program before 
attempting to solve the problem. 

For using the CHECK program prepare the following 
control cards. 

//XXXX0000 J0B (0000, 0000FT, XX00) , 'NAME' ,TIME=1 
/ /J0BLIB DD DSN=F060 9 . BAKH , DI SP = SHR, UNI T= 2314 ,V0L=SER= DUFFY 
//G0 EXEC PGM= CHECK ,REGION=100K 

/ /FT06F001 DD SYS0UT=A, DCB= (RECFM=FBA , LRECL=1 3 3 , BLKSI ZE=332 5) , 
// UNIT=SYS0UT 

//FT05F001 DD * 

where the first card is the regular FORTRAN job card used 
at this institution. 

Prepare your deck as Fig. 18 and it may be read in from 
the hot card reader. 



/* 



i 

Data deck 



Control cards of step 28 



Job card 



Figure 18. Set-up for Using CHECK Program. 
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If the output of the CHECK program has any error messages 
or has not been run, then there are some mistakes in the 
data deck or control cards. 

If there are no error messages on the output of CHECK 
program, then study carefully the output and make sure it 
is the same problem that has been intended to be solved. 

Once the correctness of the input data has been insured, 
proceed to the next step. 

Step 29 : 

AXITTS stands for Axisymmetric Transient Thermal Stress 
and this is the name given to the program presented in 
Appendix B when stored in the computer. For using that we 
must prepare the following control cards. 

//XXXX0000 J0B (0000, 000FT,XX00) 'NAME' ,TIME=10 

// J0BLIB DD DSN=F0609. BAKH ,DISP=SHR,UNIT=2314 ,V0L=SER= DUFFY 

//G0 EXEC PGM=AXITTS ,REGION=425K 

/ /FT06F001 DD SYS0UT=A , DCB= (RECFM=FBA , LRECL=133 , BLKSI ZE = 332 5) , 
// UNIT=SYS0UT , SPACE= (CYL , (6,1)) 

//FT05F001 DD * 

where the first card is the regular F0RTRAN job card. It is 
advised to ask for 10 minutes time, i.e., TIME=10, since 
this would not affect the priority of the job within the 
class K jobs. 

Prepare the deck as in Fig. 19 (it may be read into the 
computer from the hot card reader) and submit a so-called 
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service request card, available in the computer center, to 
the operator on duty. 



/* 



Data deck 



Control cards of step 29 
Regular F0RTRAN job card (TIME=10) 



Figure 19. Deck Set-up for Using AXITTS Program. 

Step 30 : 

If the user wishes to have a listing of the program he 
may prepare the following control cards as in Fig. 20 since 
the program is also listed on the data cell for this purpose. 

Regular F0RTRAN job card 
//PRNT EXEC PGM=IEBPTCH 
//SYSPRINT DD DUMMY 

//SYSUT1 DD DSN=F0609 . BAKH,DISP=0LD,UNIT=2321 , 

// V0L=SER=CELOO3,DCB= (RECFM=FB , LRECL=80 ,BLKSIZE=2000) 
//SYSUT2 DD SYS0UT=A, SPACE= (CYL , 6) 

//SYSIN DD * 

PRINT TYP0RG=P0 ,MAXFLDS=1 , M AXN AME = 1 
MEMBER NAME = AXITTS 
REC0RD FIELD= (80) 

/* 

Figure 20. Deck Set-up for Obtaining a Listing 
of the Program AXITTS. 
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The set-up deck of Fig. 20 may be read in from the hot card 
reader to get a listing of the program. 

Step 31 : 

If the user wishes to obtain a deck of the program AXITTS 
he may prepare a deck as in Fig. 21 and read it in from the 
hot card reader. He must notify the operator on duty to be 
prepared for punching almost two boxes of IBM cards. 

Regular F0RTRAN job card 
//PUNCH EXEC PGM= IEBPTCH 
//SYSPRINT DD DUMMY 

// SYSUT1 DD DSN=F- 0909. BAKH ,DISP=0LD,UNIT=2321, 

// V0L=SER=CELOO3 , DCB= (RECFM=FB, LRECL=80 , BLKSI ZE=2000) 

// SYSUT2 DD SYS0UT=B , SPACE = (CYL , 6) 

//SYSIN DD * 

PUNCH TYP0RG=P0 ,MAXFLDS=1 ,MAXNAME=1 
MEMBER NAME= AXITTS 
REC0RD FIELD= (80) 

/* 



Figure 21. Deck Set-up for Obtaining a Punched 
Deck of the Program AXITTS 

Step 32 : 

For using the program presented in Appendix B, if the 
user does not have access to the computer facility at N.P.S., 
he must punch a copy of the program. (Good luck!) 
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For the storage location 425K bytes are required and 
since the output may be longer than what usually is allowed, 
three additional cylinders are recommended. The execution 
time required depends on the size of the problem and the 
number of time integration steps; however, 10 minutes 
computer time would be sufficient for a problem of 140 nodes 
and 500 steps of time integration, with 20 temperature 
vectors stored for stress calculation. 
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APPENDIX D 



PROGRAMMING 

The computer program presented in Appendix B is formed 
from one main program and fifteen different subroutines. 

The communication between the main program and the sub- 
routines is handled through the common blocks. In order to 
minimize storage requirements, most of the storage location 
of the system stiffness matrix is used for temperature 
calculations. This is accomplished by use of EQUIVALENCE 
statements. The total storage requirement is 425K bytes. 

In this section the function of some parts of the pro- 
gram is discussed and the assumptions used are brought to 
attention. 

1. MAIN PROGRAM 

The main program is simply calling different subroutines 
when they are needed. The subroutines which are called in 
main program are: 

INPUT, PR0B , CANDY, FL0W, F0RMV , TEMPER, PL0T , STIFF, F0RMF , 
CENTF , PRESS, DISPL, AND STRESS. 

2. SUBROUTINE INPUT 

In this subroutine the data regarding the geometry of 
the body, subdivisions into elements, and the material 
properties are'read in and printed out. Calculation of the 
mid-side coordinates based on the straight line is done 
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here. This subroutine calls for subroutine PL0TRZ once 
to produce a graphical representation of the nodal points 
then the half-band-width of the system stiffness matrix is 
calculated from the connectivity array. Also, depending 
upon the choice of the system of units, the necessary con- 
versions in each system of units are accomplished and the 
reference temperature (70°F or 20°C) based on the system of 
units is selected here. 

3. SUBROUTINE PR0B 

For each problem this subroutine is called once by the 
main program to read in and print out the information 
regarding the nature of the problem. Based upon the given 
information the type of the problem is distinguished here 
and for the thermal problems the length of time integration 
is calculated. 

4. SUBROUTINE CANDY 

Subroutine CANDY (C and Y) evaluates the capacitance 
matrix C and admittance matrix Y + of Eq. 13, in banded form. 
For each problem this subroutine is called once. The func- 
tional flow chart of the subroutine CANDY is given in 

Fig. 22. Since the only non-zero elements of the symmetric 
* 

matrix Y (Eq. 11) are the diagonal, first and second super- 
diagonals, in order to conserve storage and reduce the num- 
ber of arithmetic operations, the non-zero elements are 
stored in three different vectors (DIAG, 0FF1, 0FF2) . 
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Figure 22. Functional Flow Chart of CANDY. 
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Once the system capacitance and admittance matrices C 
and Y + are formed, then the matrices A and G (Eq. 29) are 
evaluated and replace C and Y respectively. 

5. SUBR0UTINE FL0W 

If there is any convection or constant temperature 
thermal boundary condition, this subroutine is called from 
the main program for each step of time integration in order 
to evaluate the temperature of the fluid nodes. The calcu- 
lation is based upon the constant fluid flux assumption 
(for both inside and outside flow). Consider a section of 
an irregular cylindrical pipe as Fig. 23 and focus attention 
on element i+1 on the inner boundary. 




Figure 23. Fluid Node Representation for Inside Flow. 

If fluid is moving from left to right we may number the R 
and Z coordinates of the fluid nodes of the element i+1 as 
in Fig. 23, then the volume of fluid pumped into the pipe 
at time t is given by 
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Q • = a v t , 



(56) 



where a and v are, respectively, the area and velocity at 
the entrance section. Also the volume of the pipe up to 
the mid-node of element i+1 is 



where Vh is the volume of the pipe between the entrance 
section and element i+1 on the inner boundary. 

Now we can write 

3 

R = £ R. N. , (58) 

i= 1 1 1 

where are the one dimensional shape functions (given in 
Appendix A) . 

If we substitute Eq. 58 into 57 and make the coordinate 
transformation, after integration we get 



Now by equating Eq. 59 to Eq . 56 at any time t we may deter 
mine whether the front of the flow has already passed the 
mid-node of the wetted side of the element i+1. 

A similar argument can be carried out for the end node 
of the element i+1. In that case we also get an equation 
similar to Eq. 59 with different numerical coefficients. 

For the case of the outside flow it has also been 
assumed to have a constant fluid flux and a fictitious 



Q-^ - Vh + / irR dZ. 



(57) 



1 




entrance circular cross-sectional area (whose radius is the 
largest R plus unity) has been assumed. 

The numerical coefficients of R's in Hq . 59 are used 
in subroutine FL0W and it has been assumed that the temper- 
ature of a fluid "particle" does not change during passage 
through the active section. This is believed to be an 
acceptable approximation for representative values of flow 
velocity and active length. With this assumption, for. a 
given entry time- temperature relation (inside and outside 
flow) this subroutine evaluates the fluid nodal temperatures 
at any time. 

6. SUBR0UTINE F0RMV 

At every step of time integration this subroutine is 
called by the main program to form the vector v of the 
right-hand side of the discretized finite element (Eq. 13) 
for the given thermal boundary conditions. The program 
itself is self explanatory. 

7. SUBR0UTINE TEMPER 

The nodal temperatures are calculated in this subroutine 
with the trapezoidal time integrations. Irons' correction 
is applied here after every 10 steps of time integration. 

The temperature vectors selected for the stress analysis are 
stored. The transient temperatures will be printed out at 
the desired interval of time. 
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8. SUBR0UTINE STIFF 



For any stress or thermal stress problem subroutine 
STIFF is called once by the main program to. evaluate the 
stiffness matrix of the system. Once the stiffness matrix 
of the system is calculated then the desired structural 
boundary conditions are applied and, at the end, the Cholesky 
decomposition is performed. The functional flow chart of 
subroutine STIFF is given in Fig. 24. 

9 . F0RMF 

In subroutine F0RMF the thermal load vectors are calcu- 
lated for as many as (IVEC) given temperature vectors. 

These thermal load vectors are the columns of a rectangular 
matrix F. The provision is made that no matter how many 
temperature vectors are given (always IVEC £ 20) the IVEC 
number of columns of the F are filled with the thermal load 
vectors and the last column of F will contain the load 
vector corresponding to the unit end displacement for zero 
axial force when the plane-end boundary condition is applied. 
The flow chart of subroutine F0RMF is given in Fig. 25. 

10. SUBR0UTINE CENTF 

If the system is rotating about the axis of revolution, 
this subroutine is called once by the main program to eval- 
uate the centrifugal load vector. This load vector is 
always placed in the first column after the thermal load 
vectors (IVEC+1 position) in F. 
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L=1 ,NET 




Figure 24. Functional Flow Chart of STIFF. 
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Figure .25. 



Functional Flow Chart of F0RMF. 



11. SUBR0UTINE PRESS 



If there is any pressure acting on the system, the 
pressure load vector is calculated in this subroutine and 
this vector is placed in the column next to centrifugal load 
vector, if any, otherwise it will fill the position desig- 
nated for centrifugal load vector in F. 

12. SUBR0UTINE DISPL 

The displacement vectors are found in this subroutine. 
Since the stiffness matrix is already decomposed in banded 
form, then the displacement vectors one after another are 
obtained by back and forward substitution. The principle 
of superposition is applied here and finally the axial 
force, if any, is corrected for the plane-end boundary 
condition . 

13. SUBR0UTINE STRESS 

For every problem this subroutine is called by the main 
program once to evaluate the stresses at the points 
(p=±l, £= ± — - — )of each element. The transient stresses are 

/X 

calculatcd and printed for each displacement vector. For 
each problem the maximum mean stress and the maximum octa- 
hedral shearing stress, together with the corresponding 
times and locations, are printed. 

14. LIMITATIONS 

In the process of the development of the program dis- 
cussed so far, it has been intended that all of the calcula- 
tions and storage of data occur in the computer without 



169 



using any external devices such as disks or magnetic tapes 
so as to be able to solve any sizable problem in relatively 
short time. 

For this reason the following limitations (Table VI) are 
set forth which give an overall size of 450K bytes to the 
program. 

TABLE VI 

MAXIMUM VALUES FOR PROGRAM PARAMETERS 



NBAND 


Half band-width of the system stiffness matrix 


66 


NMAT 


Total number of different materials 




5 


NET 


Total number of elements 




40 


NNT 


Total number of nodes 




149 


NCF0T 


Total number of nodes in contact with 
outside fluid 




37 


NCFIT 


Total number of nodes in contact with 
inside fluid 




37 


NNRE 


Total number of nodes of right-hand end 




37 


NNLT 


Total number of nodes constrained against 
longitudinal motion 


any 


37 


NNRT 


Total number of nodes constrained against 
radial motion 


any 


37 


IVEC 


Total number of temperature vectors stored 
thermal stresses 


for 


20 


NRAMP0 


Total number of ramps for outside flow 




15 


NRAMPI 


Total number of ramps for inside flow 




15 


NPNT 


Total number of pressure nodes 




37 
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